Nonconvex QP
Adapted from: of (Floudas et al., 1999; Section 2.2), (Lasserre, 2009; Table 5.1)
Introduction
Consider the nonconvex Quadratic Program (QP) from (Floudas et al., 1999; Section 2.2) that minimizes the concave function $c^\top x - x^\top Qx / 2$ over the polyhedron obtained by intersecting the hypercube $[0, 1]^5$ with the halfspace $10x_1 + 12x_2 + 11x_3 + 7x_4 + 4x_5 \le 40$.
using LinearAlgebra
c = [42, 44, 45, 47, 47.5]
Q = 100I
using DynamicPolynomials
@polyvar x[1:5]
p = c'x - x' * Q * x / 2
using SumOfSquares
K = @set x[1] >= 0 && x[1] <= 1 &&
x[2] >= 0 && x[2] <= 1 &&
x[3] >= 0 && x[3] <= 1 &&
x[4] >= 0 && x[4] <= 1 &&
x[5] >= 0 && x[5] <= 1 &&
10x[1] + 12x[2] + 11x[3] + 7x[4] + 4x[5] <= 40Basic semialgebraic Set defined by no equality
11 inequalities
x[1] ≥ 0
1 - x[1] ≥ 0
x[2] ≥ 0
1 - x[2] ≥ 0
x[3] ≥ 0
1 - x[3] ≥ 0
x[4] ≥ 0
1 - x[4] ≥ 0
x[5] ≥ 0
1 - x[5] ≥ 0
40 - 4*x[5] - 7*x[4] - 11*x[3] - 12*x[2] - 10*x[1] ≥ 0
We will now see how to find the optimal solution using Sum of Squares Programming. We first need to pick an SDP solver, see here for a list of the available choices.
import Clarabel
solver = Clarabel.OptimizerClarabel.MOIwrapper.OptimizerA Sum-of-Squares certificate that $p \ge \alpha$ over the domain S, ensures that $\alpha$ is a lower bound to the polynomial optimization problem. The following function searches for the largest lower bound and finds zero using the dth level of the hierarchy`.
function solve(d)
model = SOSModel(solver)
@variable(model, α)
@objective(model, Max, α)
@constraint(model, c, p >= α, domain = K, maxdegree = d)
optimize!(model)
println(solution_summary(model))
return model
endsolve (generic function with 1 method)The first level of the hierarchy does not give any lower bound:
model2 = solve(2)-------------------------------------------------------------
Clarabel.jl v0.6.0 - Clever Acronym
(c) Paul Goulart
University of Oxford, 2022
-------------------------------------------------------------
problem:
variables = 33
constraints = 53
nnz(P) = 0
nnz(A) = 75
cones (total) = 3
: Zero = 1, numel = 21
: Nonnegative = 1, numel = 11
: PSDTriangle = 1, numel = 21
settings:
linear algebra: direct / qdldl, precision: Float64
max iter = 200, time limit = Inf, max step = 0.990
tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,
static reg : on, ϵ1 = 1.0e-08, ϵ2 = 4.9e-32
dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-07
iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,
max iter = 10, stop ratio = 5.0
equilibrate: on, min_scale = 1.0e-04, max_scale = 1.0e+04
max iter = 10
iter pcost dcost gap pres dres k/t μ step
---------------------------------------------------------------------------------------------
0 -1.5207e+02 -1.5207e+02 0.00e+00 5.27e-01 7.37e-02 1.00e+00 9.59e+01 ------
1 2.4388e+04 2.5071e+04 2.80e-02 7.15e-02 1.04e-02 6.91e+02 1.14e+01 9.66e-01
2 1.6512e+05 1.7080e+05 3.44e-02 3.79e-03 5.18e-04 5.68e+03 6.86e-01 9.47e-01
3 1.6216e+07 1.7985e+07 1.09e-01 7.36e-04 1.00e-04 1.77e+06 1.34e-01 8.78e-01
4 1.6199e+09 1.8067e+09 1.15e-01 1.41e-05 1.91e-06 1.87e+08 2.55e-03 9.81e-01
5 1.8390e+10 2.0786e+10 1.30e-01 2.05e-07 2.78e-08 2.40e+09 3.71e-05 9.86e-01
---------------------------------------------------------------------------------------------
Terminated with status = primal infeasible
solve time = 1.12ms
* Solver : Clarabel
* Status
Result count : 1
Termination status : INFEASIBLE
Message from the solver:
"PRIMAL_INFEASIBLE"
* Candidate solution (result #1)
Primal status : INFEASIBLE_POINT
Dual status : INFEASIBILITY_CERTIFICATE
Objective value : -1.83898e+10
Dual objective value : -2.07856e+10
* Work counters
Solve time (sec) : 1.11548e-03
Barrier iterations : 5Indeed, as the constraints have degree 1 and their multipliers are SOS so they have an even degree, with maxdegree 2 we can only use degree 0 multipliers hence constants. The terms of maximal degree in resulting sum will therefore only be in -x' * Q * x/2 hence it is not SOS whatever is the value of the multipliers. Let's try with maxdegree 3 so that the multipliers can be quadratic. This second level is now feasible and gives a lower bound of -22.
model3 = solve(3)-------------------------------------------------------------
Clarabel.jl v0.6.0 - Clever Acronym
(c) Paul Goulart
University of Oxford, 2022
-------------------------------------------------------------
problem:
variables = 253
constraints = 308
nnz(P) = 0
nnz(A) = 715
cones (total) = 13
: Zero = 1, numel = 56
: PSDTriangle = 12, numel = (21,21,21,21,...,21)
settings:
linear algebra: direct / qdldl, precision: Float64
max iter = 200, time limit = Inf, max step = 0.990
tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,
static reg : on, ϵ1 = 1.0e-08, ϵ2 = 4.9e-32
dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-07
iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,
max iter = 10, stop ratio = 5.0
equilibrate: on, min_scale = 1.0e-04, max_scale = 1.0e+04
max iter = 10
iter pcost dcost gap pres dres k/t μ step
---------------------------------------------------------------------------------------------
0 -1.3991e+01 -1.3991e+01 6.35e-16 3.32e-01 9.13e-02 1.00e+00 9.46e+00 ------
1 2.4255e+02 2.4996e+02 3.06e-02 1.45e-01 3.32e-02 8.38e+00 2.67e+00 9.00e-01
2 7.3705e+01 7.9821e+01 8.30e-02 5.20e-02 1.50e-02 6.48e+00 1.52e+00 6.49e-01
3 3.6450e+01 3.8243e+01 4.92e-02 1.15e-02 3.41e-03 1.86e+00 4.41e-01 7.97e-01
4 3.1397e+01 3.2740e+01 4.28e-02 6.37e-03 1.91e-03 1.39e+00 2.70e-01 6.34e-01
5 2.1768e+01 2.2034e+01 1.22e-02 1.23e-03 3.40e-04 2.75e-01 5.65e-02 8.06e-01
6 2.2259e+01 2.2326e+01 2.99e-03 4.09e-04 1.03e-04 6.95e-02 1.91e-02 8.17e-01
7 2.1905e+01 2.1925e+01 9.14e-04 1.28e-04 3.24e-05 2.10e-02 6.08e-03 7.38e-01
8 2.1973e+01 2.1978e+01 1.91e-04 4.09e-05 1.02e-05 4.52e-03 1.97e-03 9.90e-01
9 2.2000e+01 2.2000e+01 2.57e-06 5.82e-07 1.41e-07 6.11e-05 2.82e-05 9.90e-01
10 2.2000e+01 2.2000e+01 2.82e-08 6.39e-09 1.55e-09 6.70e-07 3.09e-07 9.89e-01
11 2.2000e+01 2.2000e+01 3.49e-10 7.90e-11 1.92e-11 8.30e-09 3.83e-09 9.88e-01
---------------------------------------------------------------------------------------------
Terminated with status = solved
solve time = 12.7ms
* Solver : Clarabel
* Status
Result count : 1
Termination status : OPTIMAL
Message from the solver:
"SOLVED"
* Candidate solution (result #1)
Primal status : FEASIBLE_POINT
Dual status : FEASIBLE_POINT
Objective value : -2.20000e+01
Dual objective value : -2.20000e+01
* Work counters
Solve time (sec) : 1.27103e-02
Barrier iterations : 11This page was generated using Literate.jl.